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Abstract 

We present a field-theoretical method to obtain consistently the equations of motion 
for small amplitude condensate perturbations in a homogeneous Bose-condensed gas 
directly in real time. It is based on linear response, and combines the Schwinger- 
Keldysh formulation of nonequilibrium quantum field theory with the Nambu- 
Gor'kov formalism of quasiparticle excitations in the condensed phase and the tad- 
pole method in quantum field theory. This method leads to causal equations of 
motion that allow to study the nonequilibrium evolution as an initial value prob- 
lem. It also allows to extract directly the Ward identities, which are a consequence 
of the underlying gauge symmetry and which in equilibrium lead to the Hugenholtz- 
Pines theorem. An explicit one-loop calculation of the equations of motion beyond 
the Hartree-Fock-Bogoliubov approximation reveals that the nonlocal, absorptive 
contributions to the self-energies corresponding to the Beliaev and Landau damp- 
ing processes are necessary to fulfill the Ward identities in or out of equilibrium. It 
is argued that a consistent implementation at low temperatures must be based on 
the loop expansion, which is shown to fulfill the Ward identities order by order in 
perturbation theory. 
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1 Introduction 



The realization of Bose-Einstein condensation (BEC) [1-4] in dilute atomic 
gases in which the atoms are confined in magnetic traps and cooled down 
to temperatures of order of fraction of microkelvins has opened a new era in 
atomic physics. The experimental and theoretical effort on the understanding 
of BEC is truly interdisciplinary [5] and the beautiful demonstration of BEC in 
atomic gases has rekindled both the experimental and theoretical interests [4]. 

The theoretical description of Bose-Einstein condensation in weakly inter- 
acting dilute gases has a long history [6-10] and has accounted for many 
experimental results [1-4]. Pioneering work on the microscopic theory of a 
Bose-condensed gas was done by Bogoliubov [11], whose original observation 
that the wave function of the excitations around the condensate (phonons) 
is a linear superposition of free particle states with equal and opposite mo- 
menta has been recently brilliantly confirmed experimentally [12]. A detailed 
and comprehensive theoretical description of the homogeneous and inhomo- 
geneous BEC is offered in Refs. [6-9,5]. 

The current experiments on trapped Bose-condensed gases can measure with 
great accuracy the dynamical aspects of the collective (quasiparticle) excita- 
tions (see Ref. [4] for a thorough discussion). The earher experiments [1-3] 
were carried at low temperatures and the measurements of the dynamics of 
the condensate revealed almost undamped oscillations with frequencies in ex- 
cellent agreement with theoretical predictions based on the time dependent 
Gross-Pitaevskii equation [13,6-9,14,15]. 

More recent experiments [3,4] probe the low energy dynamics of collective ex- 
citations at higher temperatures and show evidence for large frequency shifts 
and strong damping. There are several important processes that contribute 
to damping of collective excitations. At zero temperature a quasiparticle can 
decay into two or more quasiparticles of lower energies, a mechanism that was 
originally studied by Beliaev [16] at zero temperature and by Popov [17] at 
finite temperature in a homogeneous Bose gas. At finite temperature an im- 
portant damping mechanism in the collisionlcss regime is the Landau damp- 
ing [18-20], which has recently been studied for trapped Bose gases [15] within 
a time dependent mean-field approximation. 

The necessity for a firm theoretical understanding of the dynamical aspects of 
low-lying collective excitations in the condensed phase is underscored by the 
intense experimental effort to probe them. For a homogeneous Bose gas in the 
condensed phase a classification scheme for the different types of theoretical 
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approximations to study the excitations has been put forth in a seminal ar- 
ticle by Hohenberg and Martin [21] and more recently discussed in detail by 
Griffin [6,14,8]. A very important ingredient for a consistent description of the 
dynamics of quasiparticle excitations is the Hugenholtz-Pines theorem [22], 
which is a consequence of the underlying gauge invariance [21] and similar to 
the Goldstone theorem in that it guarantees gapless single (quasi)particle ex- 
citations if a continuous symmetry is broken with short range forces [23]. Since 
the Hugenholtz-Pines theorem is a consequence of the original gauge invari- 
ance and the Ward identities that stem from it, it is important to guarantee 
that any approximation scheme to study the dynamics of low energy excita- 
tions respects this theorem [14]. Detailed studies of the equations of motion 
for the small amplitude perturbation in the condensed phase [14,15] highlight 
the necessity for a systematic treatment compatible with the Hugenholtz- 
Pines theorem. The usual approach to obtaining the equations of motion be- 
gins with taking the expectation value of the Heisenberg equations of motion 
and invoking some type of factorization of the nonlinear correlators (see, e.g., 
Refs. [14,15] and references therein). It is at this point where approximation 
schemes to treat the nonlinear correlator may potentially conflict with the 
Hugenholtz-Pines theorem and the Ward identities. 

The focus and goals. We focus on studying the real-time evolution of 
the condensate perturbations in a weakly interacting, homogeneous Bose- 
condensed gas to begin with, but our ultimate goals are to study consistently 
the nonequilibrium evolution of collective excitations and to establish contact 
with the experimental effort. In particular, we seek to obtain the oscillation 
frequencies and damping rates of condensate perturbation both for the homo- 
geneous as well as for the inhomogeneous (trapped) cases and to study the 
nonequilibrium aspects of quasiparticle excitations in a Bose-condensed gas. 

While the study of the homogeneous Bose gas may not be directly related to 
the experimental effort, a deeper understanding of nonequilibrium phenomena 
in many-body systems is intrinsically of fundamental importance and may 
pave the way to a systematic treatment of similar effects in trapped atomic 
gases. Of particular relevance for our study, at least for the homogeneous Bose 
gas, is to assess the impact of the infrared divergences in the self-energies on the 
frequencies and damping rates of collective excitations. In a strict perturbative 
expansion in the weakly interacting homogeneous Bose gas, the quasiparticle 
self-energies feature infrared divergences, logarithmic at zero temperature and 
power law at finite temperature [6,8]. Early work by Gavoret and Nozieres [24], 
Nepomnyashchy and Nepomnyashchy [25] and Popov [26,27] points out that 
a resummation of the perturbative expansion leads to a cancellation of the 
infrared divergences in physical quantities. This result in turn suggests that 
the perturbative expansion requires a resummation to extract the quasiparticle 
frequencies and damping rates unambiguously, but the resummation scheme 
invoked should fulfill the Hugenholtz-Pines theorem. To our knowledge these 
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issues have not been so thoroughly studied and whether such divergences 
require novel treatment in the case of trapped atomic gases seems to be an 
open question. 

In this article we propose a consistent formulation of the nonequilibrium dy- 
namics by implementing a field-theoretical approach to study dynamics out 
of equilibrium. The field-theoretical approach to BEC has been clearly for- 
mulated both in the operator [28] as well as in the path integral formula- 
tion [27,29]. The functional integral formulation is particularly convenient as 
it can be generalized to study nonequilibrium phenomena by implementing the 
Schwinger-Keldysh formulation [30-35] at the level of the path integrals. Re- 
cently, Stoof [36] has advanced the Schwinger-Keldysh formulation of nonequi- 
librium field theory to study the dynamics of condensation in Bose-condensed 
gases. 

Our ultimate goal is to obtain a consistent and systematic analysis of quasi- 
particle frequencies and damping rates in both the homogeneous and inhomo- 
geneous Bose-condensed gases from a nonequilibrium approach. We begin such 
program in this article by establishing a nonequilibrium framework that allows 
to extract this information directly from the real-time evolution of small ampli- 
tude perturbations around the condensate in a homogeneous Bose-condensed 
gas. 

Brief summary of results. The main results of this article, which we con- 
sider the first in a program that is devoted to setting up the nonequilibrium 
formulation, are the summarized as follows. 

(i) We combine the Schwinger-Keldysh nonequilibrium formulation [30-35] and 
the Nambu-Gor'kov formalism [6,8,28,37] for treating the Bogoliubov quasi- 
particles in the condensed phase, along with a novel method [38] introduced 
within the context of quantum field theory, to obtain the equations of mo- 
tion for small amplitude perturbations away from equilibrium around the 
condensate directly in real time. This method leads to causal equations of 
motion and allows to study the evolution of small amplitude perturbations 
of the condensate as an initial value problem [39] which describes an exper- 
imental situation. 

(ii) The novel method [38] implemented to obtain directly the equations of mo- 
tion in real time leads to a simple derivation of the Ward identities associated 
with the underlying gauge symmetry in or out of equilibrium. In equilib- 
rium these Ward identities lead directly to the Hugenholtz-Pines theorem. 
This method is therefore seen to lead to consistent equations of motion for 
quasiparticle excitations of the condensate, that fulfill the Hugenholtz-Pines 
theorem when the homogeneous condensate is in equilibrium. We show ex- 
plicitly to one-loop order that the inclusion of absorptive contributions be- 
yond the Hartree-Fock-Bogoliubov approximation are necessary to fulfill the 
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Ward identities, which, when the condensate is in equihbrium become the 
Hugenholtz-Pines theorem. We obtain exphcitly the one-loop self-energies 
including the absorptive contributions and their spectral representations 
beyond the Hartree-Fock-Bogoliubov approximation, 
(iii) We discuss nonequilibrium aspects such as the instabilities in the quasi- 
particle spectrum when the homogeneous condensate is not in equilibrium. 
We also discuss subtleties associated with the pcrturbative expansion that 
require a consistent rearrangement of the pcrturbative series in a loop ex- 
pansion, not in an expansion in the bare coupling constant. 

The article is organized as follows. In Sec. 2 we introduce the model and the 
linear response formulation to obtain the equations of motion. In Sec. 3 we in- 
troduce the Schwinger-Keldysh formulation in the Nambu-Gor'kov formalism 
to study the nonequilibrium aspects of Bogoliubov quasiparticles. In section 4 
we introduce the tadpole method, obtain the equations of motion directly in 
real time and cast them in terms of an initial value problem. We obtain explic- 
itly the retarded self-energies up to one-loop order and obtain their spectral 
representations. In Sec. 5 we derive the generalized Ward identities, present 
an alternative derivation of the Hugenholtz-Pines theorem and confirm that 
the one-loop equations of motion obtained do fulfill this theorem and highlight 
the necessity for including the absorptive parts of the sclf-cncrgics. In Sec. 6 
we highlight relevant aspects of the pcrturbative expansions and describe the 
condensate instabilities featured in the quasiparticle spectrum when the con- 
densate is away from equilibrium. Section 7 presents our conclusions and poses 
new questions. An appendix is devoted to an alternative derivation of the Bo- 
goliubov transformation, which facilitates the Schwinger-Keldysh nonequilib- 
rium formulation. 



2 Preliminciries 

2.1 The model 

The Hamiltonian of a Bose gas in the absence of a trapping one-body potential 
is given by 



where ■^(x, t) is the Heisenberg complex scalar field representing spinless bosons 
of mass m, g — ATra/m is the strength of the pseudopotential with a > being 




(2.1) 
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the s-wave scattering length and we have set h — ks — We will recover h in 
Sec. 6, where we analyze the nature of the loop expansion. The field ip{x,t) 
and its Hermitian conjugate satisfy equal-time commutation relations 

[V'(x,i),V'^(x',i)] = 5(^)(x-x'), 
[V'(x,i),V'(x',i)] = [V'^(x,t),V'^(x',t)] = 0. (2.2) 

A more systematic treatment of interactions requires a T-matrix formula- 
tion [36] but in order to simplify the discussion we will only consider the local 
s-wave interaction at this stage. It will become clear that the method can 
be straightforwardly generalized to include the T-matrix resummation of the 
interaction. 

The Hamiltonian H is invariant under the U{1) gauge transformation 

V^t(x,t)^e-V(x,t), (2.3) 

where ^ is a constant phase. A consequence of this U{1) gauge symmetry 
is conservation of the number of particles. Indeed, the number operator of 
particles 

N ^ J d^xip\x,t)'i/j(:sc,t) (2.4) 

commutes with H and hence is a constant of motion. However, it is conve- 
nient to work in the grand-canonical ensemble in which the grand-canonical 
Hamiltonian is given by 



K = H -fxN 

+ ^ J d'xij\^,t)ij\^,t)ij{^,t)ij{^,t), (2.5) 

where the chemical potential /i is the Lagrange multiplier associated with 
conservation of number of particles. The corresponding Lagrangian (density) 
is given by 

C[i^\ ^^^+^ + f}j^p- Iv^V^^V^^ (2.6) 

where I" = ^ - ^. 

In a Bose-condensed gas, the condensate plays a crucial role [11] and hence it 
is convenient to decompose the field into the condensate and noncondensate 
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parts [14] 



V'(x,t) = (/.(: 



x,t) + x(x,t), (x(x,i)) = 0, 



(2.7) 



where 0(x, t) = (■(/'(x, t)) is referred to in the hterature as the condensate wave 
function and x(x, t) is the noncondensate operator. In the above expression, 
(0(x, t)) = Tr[pO(x, t)]/Trp denotes the expectation value of the Heisenberg 
operator C(x, t) in the initial density matrix p. The presence of the conden- 
sate 0(x, t) 7^ leads to spontaneous breaking of the U{1) gauge symmetry 
and resuhs in a profound consequence for the spectrum of the quasiparticle 
excitations. 

In the absence of trapping potentials and explicit symmetry breaking external 
sources, the condensate is homogeneous (i.e., space-time independent) and 
denoted by 0(x, t) — (po, which is the situation under consideration in this 
article. 

2.2 Real-time relaxation in linear response 

The goal in this article is to obtain directly in real time the equations of 
motion for small amplitude perturbations of the homogeneous condensate in 
an initial value problem formulation. Our strategy to study the relaxation of 
the condensate perturbation as an initial value problem begins with preparing 
a Bose-condensed gas slightly perturbed away from equilibrium by applying 
external source coupled to the field. Once the external source is switched off, 
the perturbed condensate must relax towards equilibrium. It is precisely this 
real-time evolution of the nonequilibrium condensate relaxation that we aim 
to study in this article. 

Let ry(x, t) be an external c-number source coupled to the quantum field 
'0(x, i), then the Lagrangian given by (2.6) becomes 



The presence of external source will induce a (linear) response of the system 
in the form of an induced expectation value 



where (C(x, t))^ denotes the expectation value of 0{x,t) in the presence of 
external source and 5(x, t) is the space-time dependent perturbation of the ho- 
mogeneous condensate 0o induced by the external source. The linear response 
perturbation 5(x, t) vanishes when the external source ^^(x, t) vanishes at all 
times. This is tantamount to decomposing the field into the homogeneous con- 
densate (0o)) ^ small amplitude perturbation induced by the external source 



(2.8) 



(V'(x, t))r, = 0O + (^(X, t), 



(2.9) 
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[(5(x, i)], and the noncondensate part [x(x, i)] as 

^(x, t) = <Po + 5(x, t) + x(x, t), (x(x, t))^ = 0. (2.10) 

In linear response theory 5(x, t) can be expressed in terms of the exact retarded 
Green's function of the field in the absence of external source [28,39]. An 
experimentally relevant initial value problem formulation for the real-time 
relaxation of the condensate perturbation can be obtained by considering that 
the external source is adiabatically switched on at i = — oo and switched off 
at t = 0, i.e., 

r7(x,i) =r7(x)e^*e(-i), e ^ 0+. (2.11) 

The adiabatic switching-on of the external source induces a space-time depen- 
dent condensate perturbation 5(x, t), which is prepared adiabatically by the 
external source with a given value 5(x, 0) at t — determined by 77(x). For 
t > after the external source has been switched off, the perturbed condensate 
will evolve in the absence of any external source relaxing towards equilibrium. 
Thus, the external source i]{:K,t) is necessary for preparing an initial state at 
t = and setting up an initial value problem. This method has been applied to 
study a variety of damping and relaxation phenomena in relativistic hot and 
dense plasmas [39] and damping of photons in a strong magnetic field [40]. 

Using the decomposition (2.10) and consistent with linear response by keeping 
only the linear terms in 5 and 5*, the Lagrangian C becomes (in the presence 
of external source) 

'^[x^ x] = ^o[x\ x] + Ant[x^ xi, 

with 

^o[x\x] 

Ant[x^x] 

where we have discarded the c-number (field operators independent) and sur- 
face terms. Consistent with linear response, we have only kept linear terms in 
S, S*, which are the small amplitude departure from the homogeneous conden- 
sate induced by the external source rj. We note that the Lagrangian vC[x^x] 
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(2.12) 



i d 



= X 



^' + — + H-2g\<Ponx 

d V2 



2 dt 2m 



(0oX^X^ + 0o'xx) , 



2 1 X "J-2^* 



+ 00 (/^- c/|0o|^) +V 



- g(t>oxx^x^ - g^xx^x^ - |x^x^xx + h.c. 



(2.13) 



is obviously invariant under the gauge transformations 



0*, 6\ x\ - e-^Vo, e'^W e" V, 



(2.14) 



which, as will be seen below, is at the heart of the Ward identities that will 
lead to the Hugenholtz- Pines theorem. 

Whereas in general a gauge transformation is invoked (correctly) to fix the 
condensate 0o to be real for convenience, this choice corresponds to fixing 
a particular gauge, which in turn hides the underlying gauge symmetry. In 
order to obtain the Ward identity associated with this symmetry we will keep 
a complex condensate 4>q and analyze in detail the transformation laws of the 
various contributions to the equations of motion. 



3 Nonequilibrium formulation 

3. 1 Generating functional 

The general framework to study of nonequilibrium phenomena is the Schwinger- 
Keldysh formulation [30-35], which we briefly review here in a manner leads 
immediately to a path integral formulation. For an alternative presentation 
see, e.g., Refs. [36]. 

Consider that the system is described by an initial density matrix p and a 
perturbation is switched on at a time t^, so that the total Hamiltonian for t > 
to, H{t), does not commute with the initial density matrix. The expectation 
value of a Heisenberg operator 0{t) = U^^{t,to)0{to)U{t,to) is given by 



where U{t, to) is the unitary time evolution operator in the Heisenberg picture 



with T the time-ordering symbol. If the initial density matrix p describes a 
state in thermal equilibrium at inverse temperature /3 with the unperturbed 
Hamiltonian H{t < to) = H, i.e.. 



m)) 



TYpU-\t,to)OU{t,to) 
Trp 



(3.1) 




(3.2) 



e-'''' ^U{to-iP,to), 



(3.3) 
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00 +00 



-00-i|3 

Fig. 1. The contour C in complex time plane in the Schwinger-Keldysh formulation. 
It consists of a forward branch running from t = — oo to t = +oo, a backward 
branch from t = +oo back to t = — oo, and an imaginary branch from t = — oo to 
t = —oo — ip. The sources serve to generate the real-time nonequilibrium Green's 
functions. 

then the expectation value (3.1) can be written in the form 

- ™(to-zAto)[/-Hto,t)0[/(t,to) 

TTU{to-^f3,to) ■ ^'■^> 

The numerator of this expression has the following interpretation: evolve in 
time from up to t, insert the operator O, evolve back from t to the initial 
time to and down the imaginary axis in time from to tQ — if3. The denomina- 
tor describes the evolution in imaginary time which is the familiar description 
of a thermal density matrix. We note that unlike the S'-matrix elements or 
transition amplitudes, expectation values of Heisenberg operators require evo- 
lution forward and backward in time (corresponding to the U and f/~^ on each 
side of the operator O). 

The time evolution operators have a path-integral representation [41] in terms 
of the Lagrangian, and the insertion of operators can be systematically handled 
by introducing sources coupled linearly to the field operators [43], i.e., 

C[xKx]-^^x\x]+fx + X^j- (3.5) 

The introduction of sources also allows a systematic perturbative expansion, 
since in such an expansion, powers of operators are obtained by functional 
derivatives with respect to these sources, which are set to zero after functional 
differentiation. We note that these sources j, j* introduced to generate the 
perturbative expansion in terms of functional derivatives with respect to these, 
are different from the external sources r], rj* introduced in (2.8) to generate an 
initial value problem in linear response and to displace the condensate from 
equilibrium. 

Since there are three different time evolution operators, the forward, backward 
and imaginary, we introduce three different sources for each one of these time 
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evolution operators, respectively. Taking to — > — oo, we are led to considering 
the generating functional [38,39,42] 

Z[j+,j-,f] = TrU{-oo-if3,-oo,f)U{-oo,+oo,j-)U{+oo,-oo;j+), 

(3.6) 

where U {tf, ti\j) is the time evolution operator [see (3.2)] in the presence of the 
source j and for simplicity of notation we have not displayed the complex con- 
jugate of the sources j*. The denominator in (3.4) is given by Trp = Z[0, 0, 0]. 
The generating functional Z\j^,j~,j^\ can be written as a path integral along 
the contour in (complex) time plane (see Fig. 1) 



Z[3^ , j , f] = j Vcx^Vex exp i d^x Cc [x' , X, j 



(3.7) 



where I'cX^^cX denotes the functional integration measure along the contour 
C and 



/ d'xCc[x\x,j]= / d'xC\x^+,x^,j^] - / d'xC\x^-,x-,r] 

JC J —oo J —oo 

/—oo—i/3 
d'x£[x^^,X^,f]. (3.8) 
-oo 



with d^x = J d^x dt, etc. Because of the trace and the bosonic nature 
of the operators, the path integral along the contour C requires periodic bound- 
ary conditions on the fields. The superscripts + and — refer to fields defined 
in the upper and lower branches, respectively, corresponding to forward (+) 
and backward (— ) time evolution, while the superscript f3 refers to the field 
defined in the vertical branch running down parallel to the imaginary axis. 
The negative sign in front of the action along the backward branch is a re- 
sult of the fact that backward time evolution is determined by U~^{-\-oo, — oo) 
with U the time evolution operator. The contour source j that enters in the 
contour Lagrangian Cc in (3.7) takes the values of the sources and in 
the respective branches as displayed in Fig. 1. 

Functional derivatives with respect to the sources in the forward branch give 
time-ordered products of operators, those with respect to the sources in the 
backward branch give the anti-time-ordcrcd products of operators, and those 
with respect to the sources in the imaginary branch give the usual imaginary- 
time (Matsubara) correlation functions. While the sources j"*", j~ and in- 
troduced to obtain the correlation functions via functional differentiation are 
different in the different branches, as they generate the time-ordered, anti- 
time-ordered and Matsubara correlation functions, respectively; the external 
source ?7, the homogeneous condensate 0o, and the departure from equilibrium 
5 are c-numbers and hence treated the same in all branches. 
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Writing the Lagrangian as a free and an interaction part as £ = £0 + '^int, the 
generating functional can be written as a power series expansion in the inter- 
action part, which in turn can be generated by taking functional derivatives 
with respect to the sources j, j* by identifying 



X 



± 



5 



X 



t± 



x'' - -i 



_8_ 



(3.9) 



As a result, the full generating functional along the contour C can be written 
as 



{'i 



Z\j] =exp^i / df'xCi^tfi 



'Sj*' 6j 



Zo[j], 



(3.10) 



where free field generating functional Zo[j] is given by (3.7) with the non- 
interacting Lagrangian vCo[x^x] given by (2.13). 



3.2 Green's functions 



The calculation of ZqIJ] is facilitated by introducing the Nambu-Gor'kov for- 
malism [37,28,6]. Let us introduce the (bosonic) two-component Nambu-Gor'kov 
fields 



*(x,t) 
the corresponding sources 
J(x,i) 



xH^,t) 



*t(x,t) = [xHx,t),x(x,^)l, (3.11 



j(x,t) 
/(x,t) 

and the following 2x2 matrices 



Jt(x,i) = [7-*(x,t),j(x,t)], (3.12) 





1 









1 


a+ = 







1 


, 0-3 = 


-1 



(3.13) 



in terms of which the quadratic part of the Lagrangian given by (2.13) and 
now defined on the contour C in Fig. 1 can be written as 



d 



(3.14) 
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where the time derivative is understood to be taken with respect to the contour 
C. The equation of motion for the free Nambu-Gor'kov field ^ in the presence 
of the source J reads 



d 



^(x,t) = -J(x,t). (3.15) 



The solution of this equation of motion is given by 

^_^(x, t) = - d^x'G{^ - x', t - t')J{x', f), (3.16) 
where ^(x — x', t — t') is the Green's function along this contour and satisfies 



d 

X ^(x - x',t - t') = 5c{t - t') 5^^\x - x'), (3.17) 

with Sc{t — t') the Dirac deha function along the contour dt' Sc{t — t') = 1. 
The Green's function G(x — x', t — t') has the form 

^(x - x', f - t') ^G>{x- x', t - t')ec{t - t') + G'<(x - x', t - t')Qcit' - t), 

(3.18) 

where 6c (t — t') is the step function along the contour and G^(x — x',t — t') 
obey the homogeneous equations of motion. 

The periodic boundary conditions on the fields in the path integral, a result of 
the trace over bosonic fields in (3.6), lead to the following boundary condition 
on the Green's function 

lim ^(x - x', to - t') = lim ^(x - x', to - ip - t'). (3.19) 

to — ^— oo to — >~oo 

Since along the contour to — > — oo is the earliest time and to — i(3 is therefore 
the latest time, (3.19) entails 

lim G<(x-x',to -0 = lim G>(x - x', to - - t')> V t', (3.20) 

to — >~oo to — >~oo 

which is the Kubo-Martin-Schwinger (KMS) condition for equilibrium corre- 
lation functions [43]. 

The free field generating functional Zo[J] is now obtained by writing 

*(x,i) = ^(x,t) + *j(x,i), 

*t(x,i) = ^t(x,i) + ^1^(x,t), (3.21) 
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which leads to the result 



Zo[J] = Zo[0]exp 



-I [ d^x [ d''x'j\x)G{x-x')J{x') 



(3.22) 



where and hereafter x denotes the space-time coordinates (x, t) for simplic- 
ity of notation. The source independent term Zo[0] will cancel between the 
numerator and the denominator in all expectation values in (3.1). 

Furthermore, we are interested in computing correlation functions of finite 
real times which are defined for fields in the forward (-I-) and backward (— ) 
time branches but not in the imaginary branch. For these real-time correlation 
functions the contributions to the generating functional from one source in the 
imaginary branch and another source in either the forward or backward branch 
vanish by the Riemann-Lebesgue lemma [33,38,42], since the time arguments 
are infinitely far apart along the contour. Therefore the contour integrals of 
the soTirce terms and Green's functions in the generating functional factorize 
into a term in which the sources are those either in the forward and back- 
ward branches and another term in which both sources are in the imaginary 
branch [38,42]. The latter term (with both sources in the imaginary branch) 
cancel between numerator and denominator in expectation values and the only 
remnant of the imaginary branch is through the periodic boundary conditions 
along the full contour in the Green's function. 

Thus the generating functional for real-time correlation functions simplifies to 
the following expression, defined solely along the forward and backward time 
branches [33-35,38,42], 



/ + 00 , 
d^xiC;^,[-i5/5r+, -i5/5J^ 
-oo ^ 



-/:i^t[is/sr~,is/sj-]) 



exp 



oo r+oo 



d^xd^x' 



X 



2 J — oo J — oo 

J^+(x)G'++(x,x')J+(x') + J^-{x)G—{x,x')J-{x') 



-J^+{x)G+-{x,x')J-{x') - J^-{x)G-+{x,x')J+{x') 



(3.23) 



with 



G++(a;,x') = G>{^~^,t-t')Q{t^t') + G<{^-^,t-t')Q{t' -t), 
G—{x,x') = G>{^-x',t-t')e{t' -t) + G<{^-x',t-t')e{t-t'), 
G-+{x,x') = G>{^-^',t-t'), 

G+-{x,x') = G<{x-^',t-t'), (3.24) 
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where now — oo < t,t' < +00 and the superscripts +, — correspond to the 
sources defined on the forward (+) and backward (— ) time branches, respec- 
tively. An important issue that must be highhghted at this stage, is that 
derivatives with respect to sources in the forward (+) time branch correspond 
to insertion of operators pre-multiplying the density matrix p and derivatives 
with respect to sources in the backward (— ) branch correspond to the insertion 
of operators post-multiplying the density matrix. That this is so is a conse- 
quence of the fact that the density matrix evolves in time as U {t, to)Po^~^ {t, to) 
with U {t, to) the time evolution operator. 

These four correlation functions are not independent because of the identity 



G++{x, x') + G—{x, x') - G+-{x, x') - G-+{x, x') = 0. 



(3.25) 



The diagonal elements in G~^~^{x,x') are the normal Green's functions, rep- 
resenting the propagation of single noncondensate particles, whereas the off- 
diagonal elements are the anomalous Green's functions, corresponding to the 
annihilation and creation of two noncondensate particles, respectively. 

The functions G^{x,x'), which are solutions of the homogeneous free field 
equation of motion, are simply related to the correlation functions of the free 
Nambu-Gor'kov fields Indeed, taking variational derivatives of the free 

field generating functional Zo[J] with respect to j*^ and that make up the 
source J^, one can easily show that 



G>{x,x') 
G<{x,x') 



—I 



{x{x)x\x')) {x{x)x{x')) 

(Xt(x)xt(x')) {xK^)x{x')) 

{X^ {.x')x{,x)) {x{.x')x{,x)) 

{x\x')x\x)) {x{x')xKx)) 



(3.26) 



or alternatively 

G> (x, x') = -i{^a{x)^l{x')), G< (x, x') = -i{^l{x')^a{x)), (3.27) 

where and hereafter a,h — 1,2 denote the Nambu-Gor'kov indices. The ex- 
pectation values in the expressions above are in the non-interacting thermal 
density matrix which corresponds to the quadratic part of the Lagrangian 
(Hamiltonian), i.e., the density matrix that describes free Bogoliubov quasi- 
particles in thermal equilibrium at inverse temperature (3. 

While the matrix elements G^i^{x.x') can be obtained through the usual Bo- 
goliubov transformation to the quasiparticle basis [11,6,28,8], we present in 
the Appendix an alternative derivation of these correlation functions directly 
from the spinor solutions of the homogeneous equations of motion. We find 
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that the correlation functions in (3.26) in the continuum hmit are given by 



(3.28) 



where 



G<(k,t- t')^-i[nBMgab(k)e-'^^^'-''^ 

+ [l + nB{uJk)]gab{k)e^''^'-'"^^ 



(3.29) 



In the above expressions, k = |k|, ^^(a;) = l/{e^'^ — 1) is the Bose-Einstein 
distribution function, 



m = < 



g{k) = ui 



1 -rfe0o/0o 



-rk(l>Q/4>l 



-rk(t>l/<t>o 



(3.30) 



where u^, are given by (A. 5) in the Appendix, and uJk is the energy of the 
free Bogohubov quasiparticles (see the Appendix) 



2m 



-li + 2g\(t>on -(g\(t>o?)' 



1/2 



(3.31) 



Using the relation 1 + UBiuj) = e^'^ UBiuj), one can easily verify the KMS 
condition [43] 



G> {k,t - iP - t') ^ G<{k,t - t'). 



(3.32) 



Hence, the correlation functions for the fields x^, x that will enter in the 
nonequilibrium perturbative expansion are completely determined by (3.28)- 
(3.31). 



3.3 Feynman rules 



From the generating functional of nonequilibrium Green's functions (3.23), it is 
clear that the effective interaction Lagrangian relevant for the nonequilibrium 
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calculations is given by 

^^t[x^^,X^] = Ant[x^+,x1 - Ant[x^-,X-], (3.33) 

where the fields and defined on the forward (+) and backwards 

(— ) time branches respectively. Consequently, this generating functional leads 
to the following Feynman rules that define the perturbative expansion for 
calculations of nonequilibrium expectation values. 

(i) There are two sets of interaction vertices defined by >Cf^[x^^, x^]: those in 
which the fields arc in the forward (+) branch and those in which the fields 
are in the backward (—) branch. There is a relative minus sign between 
these two types of vertices. 

(ii) There are four sets of Green's functions given by (3.24) in terms of the 
normal and anomalous correlation functions in the form displayed in (3.26). 
These correlation functions are completely determined by (3.28)-(3.31). 

(iii) The combinatoric factors are the same as in the equilibrium or imaginary- 
time (Matsubara) formulation. 



4 Releixation of condensate perturbations: An initial value problem 



The equations of motion for the small amplitude condensate perturbation 5{x) 
induced by the external source ri{x) is obtained by implementing the tadpole 
method presented in Refs. [38,42]. This method begins by writing the original 
Bose fields in the forward (+) and backward (— ) time branches as 

^^{x)^<f>o + S{x) + xH^), (4.1) 

where 0o is the homogeneous condensate in the absence of external source, 6{x) 
is the perturbation of the condensate induced by the external source which 
vanishes in the absence of external source and x^(^) the noncondensate 
fields in the forward (+) and backward (— ) time branches. The external source 
ri{x) arc c- number fields and hence taken the same value in both forward and 
backward branches. The strategy to obtain the equations of motion for small 
amplitude condensate perturbation 6{x), 6*{x) is to consider the linear, cubic 
and quartic terms in xi^)^ xK^) perturbation theory and impose the tadpole 
condition 

{xH^))r,^{x'H^))r,^0 (4.2) 

order by order in the perturbative expansion, but, consistent with linear re- 
sponse, only keep contributions linear in S{x), S*{x). 
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4-1 Lowest order (classical) equations of motion 



To illustrate this technique within the simplest setting and before embarking 
on the lengthy calculation at one-loop order, let us focus on obtaining the 

equations of motion for S{x), S*{x) to lowest order in g by fulfilling the tadpole 



condition (x^(0)). 



to lowest order in the interaction. Since the lowest 



order interaction terms in >Cint[x^x] are those linear in x^ and x [see (2.13)], 
one obtains 



/d^x|(x+(o)x^n^)) 

+ 00 [n-glM^] +v{x) 



(9 \ 



+ (x+(o)x+(^)) 



d 



- 2g\<Po\'j S*{x) - g<P*o'S{x) + 0^ (// - g\(Po\') + V*ix 

-(x+(o)x^-(x)) 



(9 \ 



+ 00 [i^ - g\(j)o\^) +ri{x) 



{X^m-{x)) 



d 



- 2^|0o|') 5\x) - gct)*^5{x) + 0S - 9\M') + | = 0, (4.3) 

which lead to the lowest order equations of motion for 5(x, t) and 5*(x, t) 



^^ + ^ + /^-2(?|0orj 5{^,t)-g4>l5*{^,t) 



+ 



00 {iJ' - g\(t>G?) +v{^,t) = 0, 



+ ^ + - 2^|0or) 5*(x, t) - g<j)*o'6{^, t) 



+ <l>l{fi-g\<Po\')+V*M^O. (4.4) 



The equations of motion obtained from the tadpole condition {x~{^))ri = 
arc the same as those given in the above expressions. This set of equations is 
recognized as the Gross-Pitaevskii equations [13,6,7] for the condensate wave 
function linearized for a small amplitude perturbation around a homogeneous 
condensate 0o- 

In order to generalize the features gleaned from the lowest order equations of 
motion to higher orders in the loop expansion, we rewrite the above equations 
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in the following illuminating form 

+ ^ + <^(x, t) - Eff 5(x, t) - Egr (x, t) + r[0o, 0^] = -r?(x, t), 

(4.5) 

where the normal and anomalous self-energies are given by 

Eff = Eg) = 2^|0or, Eg = Eg)* = (4.6) 

respectively, and the tadpole T[0o, 0o] denotes terms independent of the con- 
densate perturbation 5 and S* 

m,0o]=0o(/^-5|0o|'). (4.7) 



Before proceeding further, we note that whereas this lowest order calculation 
is fairly straightforward, the equations obtained above feature the following 
important aspects that will be shared also by the higher order calculations. 

(i) The last two lines in the tadpole condition (4.3), which are obtained by 
replacing —>■ and x"*" — > x~ the correlation functions, lead to the 
same equations of motion as those obtained from the first two hues. This 
is a generic feature to all orders in the perturbative expansion [38] and is a 
consequence of the fact that the coefficients of the correlation functions are 
local functions of t which must be the same in both branches [38,42]. 

(ii) If the homogeneous condensate |(/)oP 7^ A*/ 5' then the tadpole T[(/)o, (p^] in 
(4.4) acts as a driving force for 5, which will therefore be nonvanishing. Hence, 
the requirement that the departure from equilibrium 6 vanishes in the absence 
of external source, i.e., when 1] = at all times, implies that the homogeneous 
condensate must fulfill the lowest order equilibrium condition 

li-glM^^O. (4.8) 

In other words if and only if the homogeneous condensate fulfills the equilib- 
rium condition (4.8), the perturbation 5 will vanish when the external source 
77 vanishes at all times. 

(iii) Setting 6, r] = Oin the interaction Lagrangian (2.13) the tadpole condition 
(x^(0))j? — leads to the equation for the homogeneous condensate 

r[0o>S] = 0o(/x-^|0or)=O, (4.9) 
which determines the equilibrium value of the homogeneous condensate. 
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Now consider adding an external space-time independent source term rj, which 
leads to a space-time independent shift of the condensate (po ^ (po + S. With 
the homogeneous condensate now being 4>o + 6 the corresponding equation of 
motion is therefore given by 

T[<l>o + 5,<Pl + 5*]^-r). (4.10) 

Expanding to linear order in 6 and 6*, this equation becomes 

TM^?Ill^S+?%^S- = -,. (4.11) 

which, upon comparing to the first equation of motion in (4.5) for a space-time 
constant perturbation, leads to the identification 

dT[(j)o,(P*] _ (0) gT[0o,0*] _ (0) , . 

— — - ^11 , ^— — -^12 ■ K^-^'^) 

Furthermore, if the condition (4.9) is fulfilled, then using the expressions for 
the self-energies (4.6) we obtain 

= Eff - f Eg. (4.13) 

00 

Assuming that the condensate 0o is real, which can always be made by a gauge 
transformation (2.14), the above expression is recognized as the Hugenholtz- 
Pines relation [22]. 

In general, for an arbitrary phase of the condensate, the lowest order expression 
for the Hugenholtz-Pines theorem (4.13) accounts for the fact that E^ must be 
complex by gauge invariance. Indeed, E^^^ must be invariant under the gauge 
transformation (2.14) since it multiplies 5, whereas E^^^ must transform as an 
object with charge +2 since it multiplies S* and the full equation of motion for 
S must transform as 6 itself. Similarly, E2i^ must transform as an object with 
charge —2 under the gauge transformation (2.14). This is an important point 
that will be taken up again in higher order calculations and will be the basis 
for the alternative derivation of the Hugenholtz-Pines theorem to all orders in 
perturbation theory (see Sec. 5). 

While these features emerge as rather trivial relationships that are gleaned di- 
rectly from the lowest order (classical) equations of motion, they will manifest 
themselves to all orders in the perturbation expansion in the tadpole method. 
The much less trivial statement that follows to all orders in the perturbative 
expansion implemented via the tadpole method and that will be confirmed by 
an explicit calculation below, is that this method leads to a gapless approxi- 
mation in the classification of Hohenberg and Martin [21]. 
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(b) 



Fig. 2. Feynman diagrams contributing to the self-energies (a) Tif-^{x — x') and 
(b) E^2^(rE — x') up to one- loop order with the combinatoric factors denoted in 
parenthesis and the labels representing the time branches suppressed. A solid line 
with arrows denotes the noncondensate particle, a wiggly line denotes the conden- 
sate 4)Q or and a dashed line with open (closed) circle denotes the condensate 
perturbation 5 {5*). 

Furthermore, the lowest order equilibrium condition (4.8) relates the equilib- 
rium condensate at this order to the chemical potential by /i = 1 0o P > which is 
the lowest order Hugenholtz-Pines relation (4.13) and hence leads to a gapless 
spectrum for free Bogoliubov quasiparticle excitations in equilibrium conden- 
sate 

1/2 

(4.14) 

as can be easily verified from (3.31). 

4-2 Equations of motion at one-loop order 

The simple example above highlights the tadpole method to obtain the equa- 
tions of motion for 5{x) and 5*{x), which we now pursue to one-loop order. 
The strategy is the same, but at one-loop order we now face many Feynman 
diagrams by treating the effective Lagrangian [x x^*"] in perturbation 
theory with the Feynman rules described above. Just as in the simpler lowest 
order example worked out above, there are two types of contributions to the 
equations of motion: (i) contributions that are proportional to 5, 5*, these give 
the corresponding self-energies, and (ii) contributions that are independent of 
(5, 5*. The latter are the higher order generalization of the tadpole T[0o,0o] 
and its complex conjugate in (4.5). 

Using the Feynman rules described in the previous section, the tadpole con- 
dition (x^(0)),7 = leads to the following expression 



2m 



+ 9m\ 
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+ Ei2(x-x')5*(a;') + r[0o, 0^] + r/(a;) 



+ (x+(0)x+(a;)) 



Eii(x-x')(5(x') 



dt 



+r*[0o,0^] + 77*(x) 



d X' 



^22{X - X')5*{X') + E2l(x - 



(4.15) 



where ^^^(a; — x') = H'^j^^x — x') + S^^^(a; — x') are the self-energies and 
T[0O5 0o] ^'^^ tadpole. The diagrams for the self-energies Sj,^ and 25*2'^ up 
to one-loop order are depicted in Fig. 2 with the labels representing the time 
branches suppressed for simplicity, those for E^^" and E]'"2~ can be obtained in 
an analogous manner. The diagrams for the tadpole T[0o, 0o] depicted in 
Fig. 3. 

Introducing the space Fourier transforms for 5{x), r]{x) and T,ab{x — x') as 



J ^4(^)e^'^•^ etc., 



(4.16) 



we find the final form of the equations of motion in momentum space to be 
given by 

(d \ /■+00 I- 

X (51Jt')l + n<i>o, 4^1] + %(i) = 0, 



— « 



dt 



\ f+°° 
— + dt' 



S22(k,t-t') <5*k(0 + S2i(k,t-t') 



(4.17) 



The equations of motion obtained from the tadpole condition {x~{0))v = 
the same as those given in the above expressions. While the above expression 

has been obtained at one-loop order, it is straightforward to conclude after a 
simple diagrammatic analysis that the structure of the equations of motion 
obtained above is general and valid to all orders. 

Figure 4, which displays the diagrams for the normal self-energy En, clearly 
shows that the self-energies have a local, instantaneous contribution (the first 
two diagrams in the first bracket), which are known as the Hartree-Fock- 
Bogoliubov (HFB) contributions, and a nonlocal, retarded contribution (the 
rest of the diagrams) that will lead to absorptive parts. A simple calculation 
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Fig. 3. Feynman diagrams contributing to the tadpole T[0O)'?^o] ^° one-loop 
order. 



(2) 



(2) 



+ (2) 



Fig. 4. Feynman diagrams contributing to the normal self-energy — x'). Those 

for the anomalous self-energy Si2(x — x') can be obtained from the diagrams in 
Fig. 2(b) by attaching the labels it in the same manner as those depicted in this 
figure. 

with the real-time nonequilibrium Green's functions obtained in the previous 
section reveals that 



Sn (k, t-t') = ™ 6{t - t') + Sfi (A;, t-t')Q{t- t') 



.12 5{t-t') + i:'^^{k,t-t')Q{t-t'). 



(4.18) 



From the explicit expressions for the self-energies or by taking complex con- 
jugation of the equation of motion for 5k (t), we find the properties 



^12( O) 



S2i(k,t-t') 
E22(k,t-t') = S^l(-k,^-^')• 



(4.19) 



Furthermore, rotational and parity invariance imply that the self-energies are 
only functions of k (as can be explicitly confirmed at one- loop order). 

The HFB instantaneous parts S^^^^ and E^2™ including the lowest order con- 
tributions (4.6) are given by 



^HFB 

^11 
no = 



2g{nQ + n) 



n 



(4.20) 



where no is the density of the condensate particles, n and m are the one-loop 
normal and anomalous densities of the noncondensate particles, respectively. 
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Using the expressions for the real-time correlation functions given by (3.26) 
with (3.29) and (3.30) in the coincidence limit (i.e., f! t) and the properties 
of the Bogoliubov coefficients given by (A. 5) in the Appendix, we find 



where in obtaining the final expression of m use has been made of (A. 5) in 
the Appendix. 

As argued above, the normal self-energies En and E22 must be invariant under 
the gauge transformation (2.14), while the anomalous ones E12 and E21 must 
transform as 0o respectively. Thus using the property (4.19) it is 

proves convenient to write 



where both E/j andEo are invariant under the gauge transformation. While 
rewriting the self-energies in this manner may seem a redundant exercise, 
the main point is to highlight and make explicit their transformation laws 
under the gauge transformation. This is an important aspect that needs to be 
addressed carefully in order to extract the Ward identities, an exact result of 
the underlying gauge symmetry to be explored below. 

The gauge invariant HFB instantaneous parts E§^^ and E§^^ can be obtained 

straightforwardly from (4.20), whereas the gauge invariant nonlocal, retarded 
parts E^ and Eq can be written in terms of their spectral representation as 



The symmetric [S{k,uj)\ and antisymmetric [A(/c, u;)] spectral functions are 
even and odd functions of a;, respectively, and found to be given by (see the 
Appendix for the explicit expressions of Uq and r-q) 




(4.21) 



Eii(A;,t-t') = J:i^{k,t-t') = J^oikt-t'), 
Si2(fc, t-t') = ^Uk, t-t') = i<Polo^ ^o{k, t - t') 



(4.22) 




(4.23) 
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SB{k,uj) 



u„u 



J[r2(2-4r, + rJ)-(l-4r, + 2r 



(27r)3 " 

X [1 + nsil^q) + nB{uJp)][5{uJ -UJq- Up) + S{UJ + LJq + UJp)], 



u, 



(2T,f - + 2r,^) - (2 - 4r, + r'^)] 

X ['^B(c^g) - nB{uJp)][6{u -UJg + Up) + 5(CJ + Uq- Up)], 



AB{k,u) 



-g^no 



(27T 



■\l [rl{2 - Avq + rl) + (1 - Avq + 2rJ) + 4rgrp] 



Aiik^u) = -g no 



X [l + nB{ug) + nB{up)][5{u - Uq - Up) - 5{U + Uq + Up)], 
'^''^ uyp[rl{l-Arq + 2r^) + {2-Ar- 



q ' ^q) 



AB{k,u) = -2g no 



X [?T'b(<^9) - nB{up)][S{u - Ug + Up) - S{u + Ug - Up)], 
(fq 



Ar-pTq] 



u 



\l[2rl-2Tp{l + rl) + ^Tp 



AL{k,u)^-2g'^no 



X [1 + ns(ci;g) + nB{up)][S{u - Uq - Up) - 6{u + Uq + Up)], 

ulul [1 + 4^1 - 2rp{l + rl) + SrpVg] 
X [nB{uq) - nB{up)][S{u - Uq + Up) - 6{u + Uq - Up)], (4.24) 

where p = k + q. The subscript B denotes the Behaev damping processes in 
which one quasiparticle decays into two and the inverse process [16,17], and L 
denotes the Landau damping processes in which the quasiparticle scatters off 
another quasiparticle in the noncondensate [18-20]. The Landau damping con- 
tributions to the nonlocal, retarded parts of the self-energies arise solely from 
the thermally excited quasiparticles and hence vanish in the zero temperature 
limit. A careful comparison of the spectral representation of the self-energies 
reveals that they coincide with those obtained by Shi and Griffin [8]. 

The tadpole term up to one-loop order is given by 



^[00, 0o] = ^0 



fj, — g ino + 2n+ -^m 



(4.25) 



with n and m given by (4.21). Setting the external source rj — 0, the equihb- 
rium condition for the homogeneous condensate T[0o, 0o] — leads to 



00 



fJ'- g (no + 2n+ 

V </'o > 



0. 



(4.26) 



In equilibrium and below the critical temperature, the condensate 0o 7^ 0, 
therefore (4.26) provides one relationship between the equilibrium condensate 
and the chemical potential. The other relationship is determined by fixing the 
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total number of particles, i.e., 

n^no + n, (4.27) 

where n is the total density of the particles. Thus, in equilibrium the two 
conditions (4.26) and (4.27) completely determine the chemical potential and 
the equilibrium condensate consistently to one-loop order. 



4-3 Nonequilibrium dynamics as an initial value problem 

In an experimental situation the dynamical evolution of the small amplitude 
condensate perturbation is studied by preparing a Bose-condensed gas slightly 
perturbed away from equilibrium by adiabatically coupling to some external 
source in the infinite past. Once the source is switched-off at time t = the 
perturbed condensate relaxes towards equilibrium and the relaxation dynamics 
is studied. As discussed above, this experimental situation can be realized 
within the real-time formulation described here by taking the spatial Fourier 
transform of external source to be of the form 

77k(i)=77ke^*e(-i), e^0+. (4.28) 

The e-term serves to switch on the source adiabatically from t = —oo so as not 
to disturb the system too far from equilibrium in the process. If at t = — oo the 
system was in an equilibrium state, then the condition of equilibrium (4.26) 
ensures that for t < there is a solution of the equations of motion (4.17) of 
the form 

Sk{t) = 4(0) e^* for t < 0, (4.29) 

where 5k(0) is related to rj^ through the equations of motion for t < 0. The 
advantage of the adiabatic switching-on of the external source is that the time 
derivative of the solution (4.29) satisfies (5k(i < 0) — > as e — > 0"*". 

Introducing auxiliary quantities n^(,(A;,t — t') defined as 

^ab{k,t-t')-^^<ik,t-t'), (4.30) 

then upon using integration by parts, neglecting terms that vanish in the 
adiabatic limit e — 0+ and assuming that 0o is the equilibrium condensate 
and hence T[0o, <Po\ = 0, the equations of motion (4.17) for i > are given by 
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+ 



dt' 



^k(t)-[sf™-nf2(^,o)]r,(t) 

nfi(A;, t - t')6k{t') + nf2(A;, t - t')S\{t')] = 0, 



+ /* dt' \ug{k, t - t')d\{t') + U^,{k, t - t')dUt'] 

Jo L 



0. 



(4.31) 



The above coupled equations of motion for tlie condensate perturbation are 
now in the form of an initial value problem with initial conditions specified at 
t = and can be solved by Laplace transform. Introducing a two-component 
Nambu-Gor'kov spinor and its Laplace transform 



Ak(t) 



(4.32) 



where s is the complex Laplace variable with Res > 0, one can write the 
Laplace transformed equations of motion in a compact matrix form as 



G-\k,s)A^{s) = -\G-\k,s) - G'-^(A;,0)lAk(0). 



(4.33) 



In the above equation, G ^{k,s) is the inverse Green's function (matrix) up 
to one-loop order expressed in terms of the Laplace variable s 



G-\k,s) 



is - /cV2m + /i - Eoik, s) -{(Po/(p*o)^o{k, s) 

-((/)^/0o)So(A;, s) -is - ky2m + /i - Eoik, -s) 



with 



E^(A;, s) = ES^'^ + Eg(A;, s), Hoih s) = Eg^'^ + Eg(A;, s), 



(4.34) 



(4.35) 



where S^(A;, s) and Sq(/c,s) are the Laplace transforms of Ti^{k,t — t') and 
Eq(A;, t — t'), respectively, 



Eg(fc,5) = 

Sg(^,s) = 



dkn 



oo /cq — is 
+°° dko 



Ssik, ko) + SL{k, ko) + AB{k, k^) + ^^(/c, ko) 



ko — is 



Asik, ko) + AL{k, ko) 



(4.36) 
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In obtaining (4.34), we have made use of the properties 

which are a result of (4.19) and (4.23). The solution of (4.34) reads 



Ak(s) = - 1 - G{K s) G-\k, 0) Ak(0), 

S L J 



where 



G{k,s)^ 



1 



D{k,s) 

is + k"^ /2m — fi + T^oik, —s) 
-{(!>* /cj>o)Eo{Ks) 



X 



(4.37) 



(4.38) 



-(0oM)So(A;,s) 
-is + k"^ /2m — iJ, + T,D{k, s) 



(4.39) 



with the denominator D{k, s) given by 



D{k,s) 



IS 



k^/2m + fJ, — T,D{k, s) is + k'^ /2m — IJ> + '^oik, —s) 



+ ^Uk,s). 



(4.40) 



The real-time evolution of the condensate perturbation Ak(t) with an initial 
value Ak(0) is now obtained from the inverse Laplace transform 



(4.41) 



where the Bromwich contour B runs parallel to the imaginary axis in the com- 
plex s plane to the right of all the singularities (poles and cuts) of Ak(s) [39]. 
We note that there is no isolated pole in Ak(s) at s = since the residue 
vanishes. 

The Green's functions G{k, s) obtained above (albeit expressed in terms of the 
Laplace variable s) agree with those obtained via the Dyson's equations [16,8]. 
However, as discussed in detail in Refs. [16,8] that in equilibrium the Beliaev 
self-energies possess infrared divergences in the limit k,s ^ and hence care 
must be taken in analyzing the analytic properties of Ak(s). 

At zero temperature these infrared divergences are logarithmic, while at finite 
temperatures, the Bose-Einstein enhancement factor enhances the infrared 
divergences to a power law. At zero temperature the infrared divergences 
that appear in the Beliaev self-energies were first analyzed by Gavoret and 
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Nozieres [24] , who showed that they cancel out exactly in the final expressions 
for physical quantities to all orders in perturbation theory. Nepomnyashchy 
and Nepomnyashchy [25] as well as Popov and Seredniakov [26] have ana- 
lyzed the infrared divergences at finite temperature and also concluded that a 
resummation of the most infrared divergent diagrams leads to a finite result. 

Having set up the description of the noncquilibrium dynamics of perturbations 
away from the equilibrium condensate as an initial value problem, we will defer 
the analysis of the real-time evolution and the infrared divergences in the self- 
energies to a forthcoming article. 



5 Ward identity and Hugenholtz-Pines theorem 



5.1 General results 



A straightforward diagrammatic analysis with the Feynman rules described 
above reveals that the generic structure of the equations of motion obtained 
via the tadpole method remains the same to all orders in perturbation theory. 
When combined with the transformation properties of 0o, S and x under the 
gauge transformation (2.14) this general form of the equations of motion allows 
to derive to all orders in perturbation theory the Ward identity for the tadpole 
T[0o, 0o]- described below, this identity is a consequence of the underlying 
gauge symmetry and in equilibrium leads to an alternative derivation of the 
Hugenholtz-Pines theorem. 

First, consider the case in which S, r] — the tadpole condition (x"^(0))^ = 
leads to 



r[0o,0a = o, (5.1) 

which is the equilibrium condition for the homogeneous condensate. For a 
space-time independent shift of the condensate 4>o ^ 4>o + 6 induced by a 
space-time independent source rj, the tadpole condition now leads to 

T[cj>o + S,cf>l + S*]^-rj, (5.2) 

which expanding to linear order in 5 and 5* becomes 

+ + = (5.3) 



We now compare this equation with the first equation of motion in (4.17), 
which has the same generic structure as the full equation of motion obtained 
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to all orders in perturbation theory, for space-time independent 5 and 5*. We 
recognize 



— — — — /I — Sii(A; — 0,00 — 0), 



'0 

" =-Ei2(A; = 0,u; = 0), (5.4) 



500 

to all orders in perturbation theory. This is obviously an all-order generaliza- 
tion of the lowest order result (4.12). 

The second important ingredient and which stems from the equations of mo- 
tion (4.17) is that under a gauge transformation (2.14) the tadpole T[0O)0o] 
transforms just as 6, 4>o and r], i.e., 

r[e^Vo,e-^Vo] = e^'m>o]- (5-5) 

Taking the gauge parameter 6 infinitesimal and comparing the linear terms 
in 9 in the above expression, we find to all orders in perturbation theory the 
Ward identity for the tadpole 

^ 00, 00 = 00 ^t; 00- 5.6) 

d(po d0o 

Therefore combining (5.4) and the Ward identity (5.6), wc obtain the following 
Ward identity which is an exact relationship between the tadpole and the self- 
energies at zero frequency and momentum 

00 [/X - En (A; = 0,w = 0)] + 0^Ei2(A; = 0,w = 0) = r[0o,0o]- (5-7) 



Above the critical temperature both the tadpole and the condensate vanish, 
thus the above equation becomes a trivial identity. However, below the critical 
temperature 0o 7^ and in equilibrium T[0o, 0o] =0, the above Ward identity 
(5.7) leads to the Hugenholtz- Pines relation 

;i = En(A; = 0, a; = 0) - ^^i2{k = 0, a; = 0). (5.8) 

00 

It is customary to choose the condensate to be real by redefining its phase via 
the gauge transformation (2.14), in which case (5.8) leads to the familiar form 
of the Hugenholtz-Pines theorem. However, for a condensate with an arbitrary 
phase, the anomalous self-energy E12 must be proportional to 0q since in the 
equation of motion it multiplies d*, which transforms under gauge transfor- 
mations just as 00- This fact can be seen explicitly at both the lowest and 
one-loop order in the expressions for the respective anomalous self-energy in 
(4.6) and in (4.22), (4.23) with the spectral functions given by (4.24). However, 



30 



the product Si25* must transform just as S or 0o therefore the phase 0o/^^o 
cancels the phase of 0o ^12- In terms of the gauge invariant self-energies 
E^) and Eo the Hugenholtz-Pines theorem (5.8) becomes the more familiar 
form 

H^Eoik ^0,u} ^0) -Eo{k ^0,u} ^0). (5.9) 

We emphasize that while the Ward identity (5.7) is an exact relationship valid 
in or out of equilibrium (corresponding to when the tadpole T[0o, 0o] vanishes 
or not, respectively), the Hugenholtz-Pines theorem (5.8) or alternatively (5.9) 
is only valid in equilibrium when the tadpole vanishes. This observation will 
become important when we discuss nonequilibrium issues below. 

5.2 Confirmation at one-loop order 

We now show explicitly that the Ward identity (5.6) is fulfilled to one-loop 
order. The tadpole up to one-loop order is given by (4.25). The Hartree-Fock- 
Bogoliubov (local) contributions to the one-loop self-energies are momentum 
and frequency independent and given by (4.20), leading to the relation 

E^ - ^S"^^ ^gL + 2n- fm] , (5.10) 

which obviously do not fulfill the Ward identity (5.6) with the tadpole up 
to one loop given by (4.25). This is a known fact that the Hartree-Fock- 
Bogoliubov approximation for the one-loop self-energies violates the Hugenholtz- 
Pines theorem. 

However, using the explicit form of the nonlocal contributions for the self- 
energies given by (4.36) with the spectral functions given by (4.24), it is a 
tedious but straightforward exercise to find that (after analytic continuation 
to real frequency u; through the substitution s — > —iu; + 0+ [38,39]) 

Efi(A; = 0, u; = 0) - ^Ef2(A; = 0, u; = 0) = 2g^m, (5.11) 

which is an exact relationship up to one-loop order, obtained for arbitrary 
II and 00- This is a purely algebraic relation which we obtained without the 
use of either the dispersion relation of quasiparticles or of the lowest order 
relationship between 0o and /i implied by (4.8). Furthermore, we have checked 
explicitly that the above relation is independent of the order of limits (i.e., k — > 
is taken either before or after u ^ 0). While the individual contributions 
feature strong infrared divergences in the limit A;, a; — > at finite temperature 
(see Ref. [8] for a discussion), they cancel exactly in the combination (5.11). 
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This cancellation is a manifestation of the general results of Gavorct and 
Nozieres [24] at least for the combination that enters in the Ward identity. 

When combined with (5.10) the above relation (5.11) leads to 

Eii(A; = 0,a; = 0)-^Ei2(A; = 0,u; = 0) = ^(no + 2ri + ^m) . (5.12) 
00 V 00 / 

Collecting the above results, we find up to one-loop order that 



^[00, cf>*o] 

00 



Eii(A; = 0,a; = 0) - ^J^uik = 0,a; = 0) 

00 



H - g {no + 2n+ —m 



(5.13) 



Therefore, the Ward identity (5.6) is manifestly fulfilled up to one-loop order 
by including the nonlocal parts of the self-energies. The equilibrium condition 
T[0o,0o] = with the relation (5.12) for the one-loop self-energies shows 
that the Hugenholtz- Pines condition (5.8) is fulfilled by properly including 
the nonlocal parts of the self-energies. This in turn implies a gapless spectrum 
for quasiparticle excitations [22,14,6]. 

This is an important advantage of the tadpole method of nonequilibrium field 
theory: the equations of motion obtained at a given order in the loop expan- 
sion are causal and guaranteed to fulfill the corresponding Ward identities 
to the given order. In the classification of Hohenberg and Martin [21] the 
tadpole method leads, therefore, to a gapless approximation that fulfills the 
Hugenholtz-Pines theorem. 



6 Issues in and out of Equilibrium 



6. 1 Issues in equilibrium: nature of the perturbative expansion 



A relevant question that emerges at this stage is the nature of the perturbative 
expansion. Taking the perturbative expansion in terms of the coupling g (or 
alternatively in terms of the T-matrix [36,8]) there seems to be an inconsis- 
tency in the calculation: the HFB local contributions to the self-energies as 
well as all of the tadpole contributions are formally of order 0{g), while the 
nonlocal contributions are formally of order 0{g^). Obviously there seems to 
be a mismatch in the orders of coupling constant, but from the confirmation 
of the Ward identity at one-loop order, this mismatch is required to fulfill the 
Ward identity. 



32 



The failure of the perturbativc expansion in terms of g can also be seen as 
follows. Using the lowest order Hugenholtz-Pines relation fi — gn^ = [see 
(4.8)] so that for fixed chemical potential, according to the grand-canonical 
formulation, griQ = n ~ 0{1). The nonlocal contributions to the self-energies 
[see (4.24)] are all of order g{gno), therefore of order 0{g), hence of the same 
order as the HFB contributions and the tadpole. This mismatch of powers of 
coupling constant and the subtleties associated with this had been recognized 
in Refs. [8,44]. 

As emphasized above with the explicit calculation, the framework that leads to 
a consistent expansion, at least for temperatures well below the critical, is the 
loop expansion. The loop expansion is formally an expansion in powers of h and 
not oig. The basis for the loop expansion is that powers of h describe quantum 
corrections, which at finite temperature are mixed with thermal contributions 
and can be systematically implemented as follows. 

Firstly the Lagrangian should be formally considered to be independent of fi, 
i.e., classical. Although the coupling g depends on h through its relation with 
the s-wave scattering length, a systematic implementation of the loop expan- 
sion requires that £ ~ C(l) in terms of h and hence g must be understood as 
a parameter in the Lagrangian of order 0{1). In the path integral representa- 
tion, each path is weighted by the exponential of the action divided by h, so 
in the path integral the Lagrangian enters in the form ~ exp {j^ J d'^xCj . The 
h can be absorbed in a redefinition of the quantum fields 

X, ^ \/^X, \/^X^ (6.1) 

The path integral measure is multiplied by an overall factor which cancels out 
in all correlation functions. 

Writing the Lagrangian as in (2.12) and gathering the terms linear, quadratic, 
cubic and quartic in x, X^ iii ^^int ^-s 

r. - r(^^ + r(^^ + r(^^ + r^^^ (r o] 

then upon rescaling the fields as in (6.1) one finds 

f - >Co + + 42 + v^42 + . (6.3) 

It is now straightforward to see that (i) the lowest order term in the equation 
of motion (4.4) are of order 0{1/Vh) since they arise from the hnear term C\^l 
in the Lagrangian, (ii) the tadpole and the HFB contributions to the one-loop 
self-energies are obtained from \/hL^^^ and are therefore of order 0{\fh), and 
(iii) the nonlocal contributions are obtained by one insertion of C\J^ and one 
insertion of C^^^_ and are therefore also of order 0{\/h). Hence, the tadpole. 
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HFB and nonlocal contributions to the self-energies are all of order 0{Vh), 
therefore 0{h) as compared to the lowest order (classical) contributions. As 
a result of this counting, diagrams with n-loops arc multiplied by /l" with 
respect to the classical terms in the equations of motion. Thus when the tad- 
pole and the self-energies are computed to the same order in the h expansion 
the Ward identity is guaranteed to be fulfilled. While the identification of the 
loop expansion with a power series expansion in h is well known in the non- 
condensed phase, the presence of new interaction vertices proportional to the 
condensate 0o in the condensed phase introduces different powers of h in the 
rescaled Lagrangian. The necessity for replacing the expansion in terms of the 
coupling by a loop expansion in terms of quantum (and thermal) corrections 
was first recognized in Ref. [44], but to our knowledge the consistency of the 
loop expansion for fulfilling the Ward identity has not been highlighted in 
detail before. 



The consistency of the loop expansion still requires the resolution of another 
important issue, which has already been recognized in Ref. [44]. In equilibrium 
and for a fixed chemical potential, the equilibrium value of the condensate is 
determined by the equilibrium condition T[0o, 4>o] = 0. As argued above, the 
expression for the tadpole T[0o, 0o] ^e recognized as the 5 independent 
terms in the equations of motion in the loop expansion. Since the tadpole is 
obtained as a power expansion in h, the equilibrium value of the condensate 00 
can also be found as an expansion in h. However, in the diagonalization of the 
quadratic part of the Lagrangian Cq, the condensate 0o enters as a parameter. 
As a result, the Bogoliubov coefficients Uk, Vk and the dispersion relation of 
quasiparticles Uk depend implicitly on h through 0o (see the Appendix) . Thus 
keeping the full 0o in the free quasiparticle Lagrangian Cq and eventually 
identifying 0o with the solution of the equilibrium condition T[0o, 4>o] = in 
the loop expansion, leads to an inconsistent expansion in h. 

The remedy for this problem is to expand the full 0o in powers of H 



00 = 00 + A, (6.4) 



where (pQ ~ 0{h^) and A = hA^^^ + /t^A'^^^ + ■ ■ ■ . Thus in the free quasi- 
particle Lagrangian Cq we now replace 0o by the zeroth order term 0q, while 
the contributions that include A are now lumped together in the interaction 
Lagrangian >Cint. The new interaction Lagrangian now becomes 
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+ 77 



+ ^ + - 2^|0o|^ J 5 - ^(/'o'^* + 00 (/X - C/100 

- ^ (0oA* + + I A|2) x^X - f (200 + A^) x^X^ 
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- 2#o<^X^X - 9(f>QSx^X^ - g<t>oxx^x^ - g^xx^x^ - |x^x^xx 

+ H.C., (6.5) 
where 0o = 0q + A in all the terms that involve 0o- 

In this manner, a perturbative expansion around the equilibrium condensate 
now in terms of the free quasiparticle Lagrangian that only involves 0o 
interactions that arc systematically treated in the loop {K) expansion is con- 
sistently cast in terms of the gapless free Bogoliubov quasiparticles as the 
nonintcracting states. The equation of motion, now consistently obtained in 
the loop expansion will involve a tadpole and self-energies that are expanded 
to the same order in h, thus guaranteeing that the Ward identity is fulfilled 
and in equilibrium the quasiparticle excitations are gapless by the Hugenholtz- 
Pines theorem. A straightforward calculation reveals that to 0{fi) our calcu- 
lation is consistent by using 0q in the Bogoliubov coefficients and dispersion 
relations in the propagators. The resulting tadpole equation can be system- 
atically solved up to 0{h) and the self-energies are unambiguously of order 
0{h) with respect to the classical term. The Ward identity is fulfilled exactly 
in the manner displayed in the previous section. 

Near the critical temperature the loop expansion fails because of the presence 
of infrared divergences. Either renormalization group [44,45] or other nonper- 
turbative techniques [46,47] must be invoked for a systematic study of both 
the static and the dynamical aspects near the critical point, which of course 
fall outside the scope of this article. 



6.2 Issues out of equilibrium: condensate instabilities 



When the condensate is in equilibrium, i.e., T[0o!0o] ~ 0' solution of the 
equation of motion for small amplitude perturbations around equilibrium will 
reveal the frequencies and damping rates of the quasiparticle excitations of the 
condensate. However, the derivation of the equations of motion allows a more 
general study and in particular to address the question of dynamics away from 
equilibrium. If the tadpole T[0o, 4>o] 7^ 0, the homogeneous condensate 0o does 
not describe a situation of equilibrium. This can be seen, for example, from the 
equations of motion (4.17), which even for ?7k(t) = features a nonvanishing 
inhomogeneity T[0o, 0o]- 



35 



To understand the nonequilibrium aspects, let us first focus on the spectrum of 
quasiparticle excitations at lowest order. The excitation energies of free Bogoh- 
ubov quasiparticles uj^ given by (3.31) become imaginary for long wavelengths 
A; w when 

(//-^|0or)(//-3^|0or) <0. (6.6) 

These imaginary frequencies correspond to instabilities of the homogeneous 
Bose gas for the values of chemical potential n and condensate (po for which 
the condition (6.6) is fulfilled. These modes cannot be treated as stable quasi- 
particles and their treatment has to be modified, in particular no occupation 
number can be defined for these modes. The equations of motion in this situa- 
tion describes a nonvanishing force term proportional to T[0o, ^o] that makes 
the long wavelength condensate perturbations to evolve in time. The insta- 
bilities for long wavelength excitations reveal the growth of the condensate 
perturbations [29,48] and are qualitatively similar to spinodal decomposition 
in self- interacting field theories [49], which must be studied self-consistently 
and nonperturbatively. This is a topic that we are currently study with the 
formulation presented in this article and on which we expect to report soon. 



7 Conclusions 

In this article we have focused on establishing a program to obtain the equa- 
tions of motion for small amplitude perturbations of a homogeneous conden- 
sate based on linear response in a manner consistent with the underlying 
symmetries. We introduced a method that combines the Schwinger-Keldysh 
formulation of nonequilibrium quantum field theory, the Nambu-Gor'kov for- 
malism of the Bogoliubov quasiparticle excitations in the condensed phase 
and a novel technique, the tadpole method, used in relativistic field theory 
to obtain the equations of motion directly in real time and consistently in 
perturbation theory. This method automatically leads to causal equations of 
motion and allow to extract the Ward identities associated with the underly- 
ing gauge symmetry that are valid in or out of equilibrium. In equilibrium, as 
can be determined by a definite condition, these Ward identities lead to the 
Hugenholtz-Pines theorem, therefore this method leads to equations of motion 
that guarantee gapless quasiparticle excitations of the condensate. 

Furthermore, this method leads to a formulation of the time evolution of per- 
turbations away from the equilibrium condensate as an initial value problem, 

thereby establishing contact with potential experimental situations in which 
a particular initial state is prepared away from equilibrium. It also highlights 
that the consistent perturbative expansion is not an expansion in the coupling, 
but a loop expansion, at least at low temperature. The loop expansion can 
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be systematically implemented and is guaranteed to fulfill the Ward identi- 
ties order by order in perturbation theory. We have obtained the equations of 
motion consistently to one-loop order and showed explicitly that the inclusion 
of nonlocal (absorptive) contributions to the self-energies that are beyond the 
Hartree-Fock-Bogoliubov approximation and correspond to the Beliaev and 
Landau damping processes are necessary to the fulfillment of the Ward iden- 
tities in or out of equilibrium. 

While our focus in this article is to present the formulation that leads to 
causal equations of motion which fulfill the Ward identities and to explore 
some of their consequences, our more overarching goal is to study the real-time 
dynamics of nonequilibrium states directly in real time both in homogeneous 
and inhomogeneous (trapped) condensates. The homogeneous case studied 
here may not bear much relevance to the experimental situation, however, the 
techniques presented in this article can be extrapolated almost immediately 
to the inhomogeneous case of trapped atomic gases. Even in the presence 
of a trapping one-body potential, the underlying gauge symmetry is global, 
i.e., the gauge transformations that leave the Lagrangian invariant are space 
independent. However, an important modification will be that the condensate 
becomes localized in the center of the trapping potential and hence breaks 
translational invariance. This also implies that the tadpole T[0o(x), 4>q{^)] now 
acquires implicit space dependence through the inhomogeneous condensate 
00 (x), which in turn will require a modification of the steps that lead to the 
Ward identities. While all of these issues require a careful study, we expect 
that the method presented here will allow us to extract the relevant Ward 
identities in the inhomogeneous case as well. 

We are currently studying these aspects with particular attention to the pos- 
sible infrared divergences and their potential impact on dynamics away from 
equilibrium and instabilities of the condensates. We expect to report on many 
of these issues in a forthcoming article. 
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A Plane wave solutions and the Bogoliubov transformation 



In this Appendix we present an alternative derivation of the correlation func- 
tions for X- directly from the plane wave solutions of the homogeneous 
equations of motion (i.e., in the absence of source) for the Nambu-Gor'kov 
field given by (3.15) 



d_ 
di 



2m 



2<7|0o|' 



1*2 



^(X,t) =0. 



(A.l) 



The plane wave solution can be written in the form 



*(x,i) = 5(A;)e-'(^'=*-''-"\ S{k) 



(A.2) 



where the two-component spinor obeys 



















_94>f ik _ 











(A.3) 



where = A;^/2m — /x+2(7|0oP- This is an eigenvalue equation with a "weight" 

'il - {9\m^- The 



(73. The eigenvalues are given by — ±0;^ with Uk - 
normalization of the positive and negative energy spinors is chosen so that 



(A.4) 



where a, (3 = 1,2 (not to be confused with the components of the spinors) 
correspond to — ±0;^, respectively. That this is the correct choice stems 
from (J3 being the "weight" in the eigenvalue equation as well as comparing 
with the limit (j)o = which is the free field case. Introducing the Bogoliubov 
coefficients Uk and Vk 



Uk = I 

UkVk = 



2ujk 

*2 



9<t>f 



2uJk ' 



1/2 



Vk = -Uk 



9^f 



Vk 00 

Uk 00 



Tk = 



^k + ik) 

9\<t>o? 



(A.5) 
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we find that the positive and negative energy spinors are given by 



Uk 
Vk 



vl 

Uk 



(A.6) 



respectively. After accounting for the interpretation of negative energy solu- 
tions as antiparticles, one can therefore write the general plane wave solution 
of the homogeneous equation of motion (A.l) as 



1 



V 





Uk 








Vk 



Uk 



i(wt;t— k-x) 



(A.7) 



where V is the quantization volume. 



At the level of second quantization, one recognizes that (A.7) is the Bogoliubov 
transformation. The operators at and aj^, respectively, annihilate and create 
a Bogoliubov quasiparticle of momentum k (energy a;*;) and obey the usual 
canonical commutation relations. The correlation functions of the Nambu- 
Gor'kov fields in the density matrix that describes free Bogoliubov quasipar- 
ticles in thermal equilibrium at inverse temperature /9 are therefore found to 
be given by 



+ nB{uJk)Gab{k)e 



iujk it-t') 



^ik-(x— x') 



(Vl>t(x',t')^a(x,t)) = -E 



nB{i^k)Qab{k)e 



+ [l + nB{uk)]Qab{k)e'^'^'-''^ 



ik-(x-x') 



(A.8) 



where riBioJk) is the equilibrium distribution for Bogoliubov quasiparticles of 
momentum k 



riBitOk) = (aiLak) 
and Q{k), Q{k) are given by (3.30). 



1 



(A.9) 
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